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ABSTRACT 

We have designed a cubic spline wavelet decomposition for the Sobolev space Hq(I) where 
/ is a bounded interval. Based on a special “point- wise orthogonality of the wavelet basis 
functions, a fast Discrete Wavelet Transform (DWT) is constructed. This DWT transform 
will map discrete samples of a function to its wavelet expansion coefficients in 0(N log N) 
operations. Using this transform, we propose a collocation method for the initial value 
boundary problem of nonlinear PDE’s. Then, we test the efficiency of the DWT transform 
and apply the collocation method to solve linear and nonlinear PDE’s. 
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1 Introduction 


Wavelet approximations have attracted much attention as a potential efficient numerical 
technique for the solutions of partial differential equations [1] - [6]. Because of their advan- 
tageous properties of localizations in both space and frequency domains [8] - [10], wavelets 
seem to be a great candidate for adaptive schemes for solutions which vary dramatically 
both in space and time and develop singularities. However, in order to take advantage of the 
nice properties of wavelet approximations, we have to find an efficient way to deal with the 
nonlinearity and general boundary conditions in the PDE’s. After all, most of the problems 
of fluid dynamics and electromagnetics, which involve solutions with quite different scales, 
are governed by nonlinear PDE’s with complicated boundary conditions. Therefore, it is our 
objective here to address these issues when designing wavelet approximations and numerical 
schemes for nonlinear PDE’s. 

We will present a new wavelet collocation method designed to solve nonlinear time evo- 
lution problems. The key component in this collocation method is a so-called “Discrete 
Wavelet Transform” (DWT) which maps a solution between the physical space and the 
wavelet coefficient space. The wavelet decomposition is based on a new cubic spline wavelet 
for Hq(I) where / is a bounded interval [1 1]. In order to treat the boundary conditions an ex- 
tra boundary scaling function <j> b(x) and a boundary wavelet 0fc(x) have been used. A special 
“pointwise orthogonality” ( see(3.7)) of the wavelet functions ipj,k(x) results in O(NlogN) 
operations for the DWT transform where N is the total number of unknowns. Therefore, the 
nonlinear term in the PDE can be easily treated in the physical space, and the derivatives of 
those nonlinear terms then computed in the wavelet space. As a result, collocation methods 
will provide the flexibility of handling nonlinearity (and, also the implementation of various 
boundary conditions) which usually are not shared by Galerkin type wavelet methods and 
finite element methods. 

The rest of this paper is divided into the following five sections. In section 2, we introduce 
the cubic scaling functions <f>(x), <f>b(x) and their wavelet functions ip(x),rpb{x). A multires- 
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olution analysis (MRA) and its corresponding wavelet decomposition of the Sobolev space 
Hq(I) are constructed using 4>b{x) and vH 1 )) ^6( x )- Then, we show how to construct 

a wavelet approximation for function in Sobolev space H 2 (I) which the solutions of PDE’s 
will belong to. In Section 3, we discuss the fast discrete wavelet transform (DWT) between 
functions and their wavelet coefficients. In Section 4, we discuss the derivative matrix T> asso- 
ciated with wavelet interpolations. In Section 5, we present the wavelet collocation methods 
for nonlinear time evolution PDE’s. In Section 6, we give the CPU time performance of the 
DWT transforms and the numerical results of the wavelet collocation methods for linear and 
nonlinear PDE’s, and a conclusion is given in Section 7. 

2 Scaling functions <^(a;), <^(z) and wavelet functions 

tp(x),tl> b (x) 

Let / denote any finite interval, say / = [0,£] and L is a positive integer (for the sake of 
simplicity, we assume that L > 4), and H 2 {I) and Hq(I) denote the following two Sobolev 
spaces with finite L 2 norm for up to the second derivatives, i.e. 

H\l) = {/(*),* € 1\ ||/"' l || ! < oo, < = 0,1,2} (2.1) 

HoU) = {/« € « ! (/)| /(0) = f(a) = SW = f'W = 0}. (2.2) 

It can be easily checked [7] that Hq(I) is a Hilbert space with the inner product 

<f,9>=J i f"(x)g"(x)dx, (2.3) 

thus, 

111/111 = \Z</./> (2.4) 

provides a norm for 

In order to generate a multiresolution for Sobolev space Hq(I), we consider two scaling 
functions, an interior scaling function <})(x) and a boundary scaling function <f>b{x) (see Figure 
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<f>(x) = n 4 {x) = 1 (?) - i)+ 

(2.5) 


M x ) = 'l x l ~ j£ x + + - l )+ - ~ 2 )+ 

(2.6) 

where N 4 (x) 

is the 4th order B-spline [13] and for any real number n 



x n f *" ^ * > 0 

+ 1 0 otherwise. 


In a pair 

they satisfy the following two-scale relationship, 


Lemma 1 

4{ x ) = H2 3 (j)<£(2x - k) 

fc= 0 



M x ) = fl-i<f>b{2x) + - k) 

L A 

(2.7) 


K — U 

q ^17 13 

here /?_i = /3 0 = --,A = » fa = — 


We summarize some properties of <f>(x) and <f>b{x) in the following lemma. 


Lemma 2 Let <j>{x) and <j>b{x) be defined as in (2.5) and (2.6), then we have 



(1) supp{<f>{x)) = [0,4]; 

(2.8) 


(2) supp((f>b(x)) = [0,3]; 

(2.9) 


(3) <*(x),<M*)€// 0 2 (/); 

(2.10) 


(4) *'(1) = -^'(3) = i, m = 0, tf(l) = 1,^(2) = -i; 

(2.11) 


(5) ^(i) = m = i m = |, hw = 

(2.12) 


For any j , k € Z, we define 


<£,,*(*) = <£(2 j £ - fc), <foj(x) = <M2 j £). (2.13) 

And for each j , let V) be the closure under norm |||/||| in (2.4) of the linear span of 
{^ iifc (a;),0 < k < 2 j L - 4,<f> b j(x),<f) bi: (L - *)}, namely 


Vj = span{<j>j tk (x), 0 < k < 2 PL - 4, 4bj(x), <f>b,j{L — x)}. (2-14) 
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Theorem 1 Let Vj,j E Z + be the linear span of (2.14), then Vj forms a multiresolution 
analysis (MRA) for Hq(I) equipped with norm (2.4) in the following sense 

(i) Vo C C Vj C • ■ • ; 

(ii) cloSffi (UjgZ+ Vj) = //£(/); 

(hi) n ; ez+ Vj = V 0 ; and 

(iv) for each j, {<f>j y k( x ) = 4>(2 3 x — fc), <)>bj( x ) = <f>b(2 Jx )i <f>bj(L — %)} is an unconditional 
basis of Vj. 

Proof. The proof for (iii) and (iv) is straightforward and omitted here. The proof for 
(i) follows from (2.7) in Lemma L In order to prove (ii), we recall a familiar result on 
interpolation cubic spline approximation for smooth functions taken from [14] and rewritten 
for the proof of our theorem, 

□ 

Lemma 3 Let 7r be the partition given by X{ = ih , 0 < i < n, h = and ,s(;r) be the cubic 
spline interpolating /(z) G (7 4 [«,6] at all points in 7r, 

■^(^t) — »/* ( ^* z) 7 

and satisfying the following boundary conditions: 

s'(a) = /'(a), s'(6) = /'(&). (2.15) 

Then s(z) uniquely exists and 

” lM r) -/ (r) ||^ <er||/ ( 4 ) IU^ 4 - r , r = 0,l,2,3 (2.16) 

where c 0 = gfj.ci = ~,c 2 = |,c 3 = 1. 

Proof of (ii) of Theorem 1 Let h — a = 0 and b = L. Consider f(x) E Cq > ( 0, L), Since 
C^°(0, L) C C 4 [0, L] H Hq( 0, L), by Lemma 3, there is an unique cubic spline corresponding 
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the partition 7 r interpolating f(x). From the fact that /( 0) = f(L ) — /'( 0) — f'{L) — 0, we 
have s(x) in Vj and then 

V - 4 

s(x) = C-i(j) biJ (x) + ^2 C k<f>j,k( X ) + c L-3<f>b,j(L - x) ( 2 - 17 ) 

k=0 

such that 

s{xi) = f{xi), 0<i<2 j L (2.18) 

where V = 2 ; L, x, = ^ j. 

Finally, from (2.16) in Lemma 3 with r = 2 we have 

|IN-/||| = lk (2) -/ {2) ll<^||/ (4) ||/2 2i . 

Therefore, as j — > oo, |||s - /||| — ♦ 0. This proves that C£°(0, L) C clos H s(U je z+ Vi). 
Then, Theorem 1 (ii) follows from the fact that C^°( 0, L ) is dense in Hq{ 0, L). 

□ 

To construct a wavelet decomposition of Sobolev space Hq(I) under the inner product 
(2.3), we consider the following two wavelet functions </>( x ),</>6( x ) ( see Figure 2) , 

t l>(x) = -^<f>(2x) + —<f>(2x — l) — -(f>{2x-2)eVi (2.19) 

th( x ) = y^<M2z) - ^<^(2x) e V\. (2.20) 

It can be verified that 

rp(n) = V»6 (n) = 0, for all n € Z. (2.21) 

Equation (2.21) will be important in the construction of the fast DWT transform later. And, 
equations (2.19) and (2.20) imply that t/>(x) and V’t(^) both belong to V\. As usual, we define 
the dilation and translation of these two functions 

x h,k{’ x ) = 4>{2 j x - k), j > 0, k = 0, • ■ • ,nj — 3, (2.22) 
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(2.23) 


Ki(*) = *(**)■ = A(2 S (L - i)) 

where rij — 2 3 L. For the sake of simplicity, we will adopt the following notations 

(■*" ) = ^b,ji x )i -2 = '4 > b,j( X )‘ (2.24) 

So when A: = — l,ny — 2, ipj t k(x) will denote the two boundary wavelet functions, not the 
usual translation and dilation of 0(x). 

Finally, for each j > 0, we define 

Wj = clos H 2 < k = -1, • • • ,nj — 2 > . (2.25) 

Theorem 2 The > 0 defined in (2.25) is the orthogonal compliment of Vj in Vy +1 

under the inner product (2.3), i.e. 

(1) Yj+i — Y> © Wj ; for j £ Z + . Here 0 stands for VfcJLWj under the inner product (2.3) 
and Vj+i = Vj + Wj. Therefore, 

(2) W J ±W j+u jeZ+; 

(3) // 0 2 (/) = W,© j6 z + W,. 

Proof. . ( 1 ) We only have to prove Vj © Wj for j = 0, namely, for 0 < / < L - 4, 0 < k < L - 3, 


< <f>( x - 0> ^i x -k)> - 0 (2.26) 

< <f>(x - > = 0 (2.27) 

< 4> b {x),il)(x - k) > = 0 (2.28) 

< f>b{x), 4>b{x) > = 0. (2.29) 


Integrating by parts twice in (2.26) and using the fact that ij>(x),4(x) £ H£(I) , we have 
< <f>(x - /), xj>(x - k) > = f <j)"(x - l)tp"(x - k ) dx 

J o 

= 4>'{ x - W(x - k ) lo - [ (jf 3 \x - l)xf\x - k) dx 

J 0 

= - / <t> [3 \x - lU'(x - k ) dx 

Jo 

= -<f> {3) {x - /)^(x - fc)|£ + / <j) (4 \x - l)xl>{x - k) dx 

J o 

= / <^ 4 l(x — /)i/>(x — fc) dx. 

JO 
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Using equation (2.21) and the identity 

4 j = o v ' 

where £(i) is the Dirac.-<5 function, so we have 

< <j>{x - l),xl>{x - k) >= j £ (j) (“I YHJ ~ ( k ~ 0) = °- 

4 3=0 V ' 

Equations (2.27) - (2.29) can be shown similarly to be true . So (1) follows from (2.19) 
and (2.20) and the fact that dimVj — 2 'PL — 3 and dimWj = 2 3 L and dimVj+i = 2 J+1 L — 3 = 
(2 J L — 3) + 2 J L = dimVj + dimWj ; 

(2) follows from (1); 

(3) follows directly from Theorem 1 (ii). 

□ 

As a consequence of Theorem 2, any function f(x ) € Hq(I) can be approximated as 
closely as possible by a function fj(x) € Vj = Vo ® Wq © W\ • ■ • ® Wj for a sufficiently large 
j, and fj(x) has an unique orthogonal decomposition 

fji x ) = fo + 9o + g\ + 1 ‘ ' + 9j (2.30) 

where / 0 € Vo, < 7 , 6 W u 0 < i < j. 
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Approximation for function in H 2 (I) 

Consider the following two functions 

t?,(x) = 2x + - 3x 2 + + ^x+ - ^(x - 1)+ (2.31) 

V2(x) = (2.32) 

For any function /(x) € // 2 (/), by the Sobolev embedding theorem we have /(x) £ C*(I) 
and, therefore, we can define the following boundary interpolation If,j/(x), j > 0 


= + Qf + ot3Tji(2 J (L — x)) + a A y 2 {y{L - x)) (2.33) 

such that 

l6j/(of= /(0), I ,J(L) = f(L) (2.34) 

l6„/'(0) = /'( 0), I b J'{L) = f(L). (2.35) 

It can be easily verified that, in order to have hjf satisfy conditions (2.34) - (2.35), we 
have to take 

®1 = - |/(°). = /(0) (2.36) 

«3 = -^-j/(i), <* = /(!). 


In many situations we do not have the values of derivatives /'(0), f'(L). However, they 
can be approximated by finite differences using only the values of /(x). To preserve the 
correct order of accuracy for a cubic spline approximation, we suggest using the following 
approximations 

/'(0) = i + <?(/,.) (2.37) 

n k= 0 

/'(£) = -T ’tckf(L-kh) + 0(V). 

n k= 0 

where h > 0 and p > 3. For p = 3, if we take 

11 

Co = - — , ci = 3 
6 
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C2 = " 2 ’ 


C3 3' 


then 5 = 3 in (2.37), and thus, equation (2.35) is satisfied within an error of 0(h 3 ). Corre- 
spondingly, the coefficients a*,, 1 < k < 4 for Ibjf(x) become 


a» = E4/(% «a = /(0) 

Jfc= 0 

cv 3 = -E4/(L-a), a 4 = /(L) 

Jfe=0 


(2.38) 


where 


, 1 3 

C ° ^‘2i +1 /i C ° 2 ’ 


Ci. = 


Cfc 


' k 2J+ 1 /fc’ 


1 < < p. 


Now we have /(x) — I(> , j / ( I ) G Hq(I) and the decomposition (2.30) can be applied for 
it. Therefore, we can find an approximation /j(x) for any function f(x) G H 2 (I) as close as 
possible, provided that j is large enough, in the form of 


fj{x) = I b,jf + fo + 9o + 9\ b 9] 


(2.39) 


where f 0 (x) G Vo, gi € 0 < i < j. 

3 Discrete Wavelet Transform (DWT) 

In this section, we will introduce a fast Discrete Wavelet Transform (DWT) which maps 
discrete sample values of a function to its wavelet interpolant expansions. Such expansion 
with the wavelet decomposition will enable us to compute an approximation of the derivatives 
of the function. 

Interpolant Operator Ivo in Vo 

Consider any function f(x) G H$(I) and denote the interior knots for V 0 by 

x[~' ) = k, fc = 1, • • • , L — 1 (3.1) 

and the values of f(x) on by 

fl~ X) = /( 4 _1) ) 5 k = !>•••, L~ 1 - ( 3 - 2 ) 
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The cubic interpolant Iy 0 f{x) of data {fj. can be expressed as follows, 

L—4 

Ivof(z) = c^<f> k (x) -f J2 Ck<f>oA x ) + - X ) (3.3) 

k - 0 


and lyof(x) interpolates data fl l \ k = !,*• 

• , L — 1, namely 


HH 

o 

i 

ii 

1 

II 

(3.4) 


Let B be the transform matrix between C ^ = (/j • • • ,/|_V) T and the coefficient 

c = (c_ l7 - • • ,c l _ 3 ) t , i.e. 

f M) = Be (3.5) 



where rij = DimWj = 2 3 L. 

It can be easily checked that the interpolation points for Vo in (3.1) and {zj^} 

for Wj,j > 0 in (3.6) satisfy a “point-wise orthogonality” condition. 

Point Orthogonality of {x^} for j > z, — 1 < k < rij — 2, 

^jA x k ] ) = 1 

0 iifc (*|' ,) ) = 0, -1 < e < n, - 2 if i > 0; 1 < i < L - 1 if i = -I. (3.7) 
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This orthogonality condition will be crucial in obtaining a fast Discrete Wavelet Trans- 
form (DWT). 

The interpolation I u , J /(x) of a function f(x) € //q(/) in Wj,j > 0 can be expressed as a 
linear combination of & = — 1, • ■ • , nj — 2, namely 

I .,/(*) = ’£ fiM*) < 3 ' 8 ) 

and 

Iu. J /(^' ) ) = /(x«), -l<k<n>-2. 

If we denote M, as the n,— th order matrix which relates f (j> = - 1 , • ■ • , and 

fh) = (/(x!fj),---,/(x«_ 2 )) T , then 

f U) = (3.9) 


where 



_L 

' 14 
1 

’ 14 


‘ 14 

1 


V 


14 


1 _J- 

14 1 


-n 1 


i_ 

14 


’ 14 


\ 
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l / 


The solution of the coefficients f hk , -1 < k < nj -2 again involves solving tridiagonal system 


(3.9) which costs (5 n-j) operations. 


Now let us assume that the values of a function f(x) € //£(/) are given on all the inter- 
polation points {x[ :) } defined in (3.1) and (3.6), we intend to find the wavelet interpolation 


Vjf{x ) € Vo © Wo ® W x • • • © Wj for J > 0, i.e. 

L — 4 

Vjf(x) = + + 

k = 0 

j=0 fc=-l 

= /-i(*) + S/j(*) ( 3,1 °) 

J=0 
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where 


7lj -2 

/-.(■) = W(x) e Vo, /i(x) = Y A«M*) € W'i.i > 0, 

k =- 1 

and 

TVO'") = /(*i' ,) ), \<k<L-\ 

Vjl(xf) = /(4 J) ), j > 0, -1 < * < n, - 2. (3.1 1) 

Let us denote f = (f( -1 ), C°\ ■ • • , f^) T the values of f(x) on all interpolation points, he. 

f'-" = {/(4“ ,, )}fci\ 
f (,) = {/(zPnET-i, J> 0 , 

and f = (f^ _1 \f*°\ • • • ,f^) T the wavelet coefficients in the expansion (3.10) 

f'-' 1 = 

f® = (An rr.„ j>o. 

The following algorithm provides a recursive way to compute all the wavelet coefficients 
f, and also the wavelet expansion (3.10) can be expanded as needed to include higher level 
wavelet spaces Wj, J + i <j< J' by adding only terms from the higher wavelet spaces, i.e. 
Wj +u ---,Wj>. 

DWT transform 
f — v f 

This direction of transform is straightforward by evaluating the expansion (3.10) at all 
the collocation points {xj^}, j > — 1 to obtain f. The “Point-wise Orthogonality” (3.7) 
of the interpolation points and the compactness of suppxf >,A:(x) can be used to reduce the 
number of evaluations. 

Number of Operations 

Let N be the total number of collocation points and N = ( L — 1 ) + f2j-o n j = 2‘ /+1 L—\. In 
the evaluation of Vjf(x^), values of and <f>(x) at dyadic points ^-,0 < k < V L,j > 0 
are needed and they can be computed once for all for future use. 
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Recalling (3.10) and the “orthogonality condition” (3.7) of the interpolation points, we 
have 

Vjf(z i _1) ) = /-i(4 _,) ). \<k<L-\ 

which needs 7 (L — 1) (flops). 

For each 0 < j < J, to compute Vjf(x^), — 1 < k < nj — 2, it needs (5j 4- 12) *n,j (flops). 
Thus, it takes 7 (L - 1) + £/ =0 (5j + 12 )n,- = 2' 7+, L(5J + 7) + 5L - 7 < 7NlogN(//ops) to 
compute the vector f. 


f — * 


f 


Recalling that f = (f (-1) , f(°), • • • , we proceed to construction of Vjf( x) in the 

following steps. 

Step 1 

Define 


f-i{x) = Ivof (_1) = + £ /-!,*&(*) + f-uL-sML - *). 

fc =0 

so /_i(x) interpolates /(x) at the interpolation points x£ — 1 < k < L — 1, namely 




(3.12) 


Step 2 
Define 

no — 2 

AM = i»o(f ,0) - (Ivo/) 10 ’) = E (- 3 - 13 ) 

/=-l 

where (W)< 0) = {Ivo/t*?’)}^, 

As a result of the “point-wise orthogonality” conditions (3.7) of the interpolation points, 
we have ipo,i( x k = 0, — 1 < / < no — 2, 1 < A; < L — 1, thus 

/o(4 _1) ) = 0, 1 <k<L-L 
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So we have 


/-i(4‘") + m 4~'’) = /~,(4~”) = w(4" n ) = /(4' n ) i <*</--! 
/-■(xfl+ZoiiD = W(4 0) ) + (<1 0) - (W)£”) = ff = /(4 0) ). (an) 


Equation (3.14) implies that function f-i(x) + fo(x) actually interpolates f(x) on both 
interpolation points for Vo and the interpolation points {x[°^} for Wq. 


Step 3 

Generally, we define for \<i<j 


/,-(*) = (3i5) 

= ”£ AmM*)- (3.16) 

Jt=- 1 

where (Vj-if){ }) = Vj-if{x[ j) ), -1 < k < rij - 2. 

Again, as in step 2 we can verify that function f-i(x) + fo(x) + ■ ■ ■ + fj{x) interpolates 
function f(x) on all interpolation points ■ {i^}. Especially, for j = J we have 

Vjf(x) = /_j(x)-(- /o(z)H 1 - fj{ x )i which will satisfy the required interpolation condition 

(3.11). 

Number of Operations. 

For j = —1, the number of operations to invert (3.9) using Thomas algorithm to obtain 
{/1 -1) } is 5 L(flops). For 0 < j < J, the cost of computing the coefficients ffp in fj(x) = 
I, Uj (/ (j) - (Vj -1 f) U) ) = £-i <k< n) --2 consists of three parts: (1) evaluation of 

— {'Pj-if( x {^)} — (^j + 7)n J (flops)-, (2) calculating the difference /0) — 

- rij(flops ); (3) inverting the matrix Mj in (3.9) - 5 rij(flops), totaling (5 j + 13 )n J {flops). 
So the total cost of finding f = 5L + £j_ 0 (5_; + 13)nj = 52j =0 (5j + 1 3)2 J < 6NlogN where 

again N = 2 J+1 L — 1. 

Now let us go back to (3.10) to see the meaning of the wavelet coefficient {fj,k} in the 
finite wavelet decomposition of space Hq(I) for function f{x). For this purpose, we first take 
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a look at the wavelet coefficient in the finite wavelet decomposition of space L 2 (/), i.e. in 
the decomposition Vj = Vo ® Wo 0 W\ • • • 0 Wj. The orthogonality here is in the sense of L 
norm not of H'q(I) norm. For simplicity , we still use the notation <f>{x) and to denote 
the scaling function and the wavelet function for this decomposition, while keeping in mind 
that they have different definitions from those of tf 0 2 (/) . Then, we can write the wavelet 
coefficient fj k in the finite Recomposition as, 

fjk = J i f(x)il>j k (x)dx 

where {V^} is the dual wavelet basis of {i>jk} in Wj. i.e.{x^* k } is such a basis of Wj that 


J *l>jk( x )M x ) 


dx = Sik 


where £/* is the Kronecker symbol. 

Using a similar method in [12], we can prove that 




where satisfies the estimate 


;=i 


i4?i < a'a''-*' 


(3.17) 


(3.18) 


with 0 < A < 1 and K a constant. 

In order to estimate f /t/>, we quote the following theorem from Meyer’s book [9]. 
Theorem A Let g(x) be compactly supported, n times continuously differentiable and have 
n + 1 vanishing moments: 



for 0 < p < n. 


Let a,0 < cv < n, be a real number that is not an integer and f(x ) € L 2 . Then f(x) is 
uniformly Lipschitz of order o over a finite interval [a, 6] if and only if for any k G Z and 
Z such that 2 ~^k € ( a, 6 ), 


| / f(x)g( Vx - k) dx | = Op-V+’W) as j — > oo. (3.19) 
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From Theorem A, we can claim that the absolute value of the wavelet coefficient 
depends upon the local regularity of f(x) in the neighborhood of the abscissa 2 ~’k. More 
precisely, if 2 ~*k E (a,£>), the decay of \fjk\ depends upon the Lipschitz regularity of f(x) 
over the interval [a, 6], as the resolution 2 J increases. This property of the wavelet coeffi- 
cients allow us to detect the location of the singularity of the function and, then, provide a 
general knowledge of the distribution of wavelet basis functions whose coefficients are larger 
in magnitude than a given threshold. The detail can be referred to [11], In the framework 
of L 2 , the wavelet function 0(x) has at least 1 vanishing moment. Hence the property of the 
wavelet coefficients mentioned above is always valid. 

Now we return to the wavelet coefficients {fjk} in (3.10). It can be easily checked that 

= ^(2 - cosw)(^^) 4 e~^r\ 

Hence 0(0) = f , which implies that 0 has no vanishing moment at all (see Figure 3) . Since 
the wavelet decomposition we considered here is in the space Hq(I), therefore, the decay 
property for the wavelet coefficients fjk ought to be related to the vanishing moments of the 
second derivative of 0(x) (seee Figure 4), not to those of 0(x). We shall illustrate this more 
precisely. 

Let {0* fc } be the dual basis of {ifjk} in Wj (recalling that the space we consider here is 
Hq) . It can be proven that for 'I'**., (3.17) and (3.18) still hold. Then we have 

A = / /"(*)«■;*)"(*) <&. (3.20) 

Notice that spline wavelets 0(x) and 0f,(x) defined by (2.19) and (2.20) are continuous 
and their second derivatives have 2 vanishing moments. Then applying Theorem A to the 
second derivatives of 0(x) and 0t(x) (namely, by taking g(x) in Theorem A to be 0"(x) and 
?/>"(x) respectively), we can prove the following. 

Lemma 4 Let 0 < cv < 1 and / E If the second derivative of the function / is Holder 

continuous with exponent cv, at Xo E /, i.e. 

|/"(x) - f"{x 0 )| < (7|x - xo|°, x E (x 0 - 8,x 0 + 8) C / 
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for some S > 0, then for any k € Z,j 6 Z + such that 2 ~*k € (xq — S/2, xq + S/2), 


\fj,k\ = 0(2 (a+1)j ), as;->oo. 


(3.21) 


Proof. We have , since ip and ip b defined by (2.19) and (2.20) are continuous and their 
second derivatives have 2 vanishing moments, by Theorem A, 

l/jtl = I </,<&>!< E l««ll</.V’,i>l+ E l««l I </.*<> 

2 “ i l£(xo— SyXQ +5) 2“ J / — 

< £ A'Al fc - ,| 0(2-( ft+1)j ) + 53 AW l '- fc| (7 = 0(2- ( “ +1)i ) 

2-]le(x 0 -6,x 0 +S) |/-A:|>2>§ 

where C in the first but last equation is a constant which depends on the second derivative 
of /(x). 

□ 


Lemma 4 implies that the wavelet coefficients fjk, j > 0, still reflect the singularity of the 
function to be approximated. In practice, when we solve PDE’s using collocation methods, 
we often use the values of the functions, not their derivatives. Therefore, in order to use the 
wavelet coefficients to adjust the choice of wavelet basis functions, we have to establish a 
relation between the magnitude of the wavelet coefficients ^ 0 and f(x). Let us first 

state the following result on the inverse of tridiagonal matrix from [15]. 


Lemma 5 Let Abeanxn tridiagonal matrix with elements a 2 , «3? * * * > a n on the subdiagonal, 
& 1? 6-2, ■ - ■ , bn on the diagonal and c 2 , c 3 , • • ■ , c n on the superdiagonal, where a,, c x ^ 0. Define 
the two sequence {u m }, {v m } as follows: 


u 0 


— 0, U\ — 1 , U in — — 2 d“ ^7?i — l^m- 1) m ^ 


U n +i — 0, V n — 1 , V m — (frm-fl ^m+1 “H C?7i+2^m+l) ^ ^ 1 

where and c n+t are arbitrary nonzero constants. Then A~ x = (cvij) is given by 


a i,j — 


FT Zk. 
aiv Q a k 

U J V ' TT £}l 
aiVQ a k 


i < j 
i > j 


(3.22) 

(3.23) 


(3.24) 
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Corollary. Let Mj be the interpolation matrix in (3.9), then we have the following 
estimates on Afj" 1 = 


!°ijl — 


K 

oM-'l 



where K = 1.1726 and a = 7 + \/l92 = 13.928. 
We delay the proof of (3.25) to the Appendix. 


Theorem 3 Let f(x) e //q(0, L) and A/ = maxj\f(x)\ and I Wj f(x) be its interpolation in 
W } defined in (3.8) and if for e > 0, — 1 < k\ < k 2 < rij — 2 


|/(*f>)| < f for ki < k < k 2 . 


then define 




where / = /(e) = mm(^, We have 




(3.26) 


(3.27) 


where C(Af) = |^(o + A/), K = 1.1726 and a = 7+ \/\92 = 13.928. 

Proof. From (3.9), we have 

f (j) = Mj’fW 

where fW = , ■ • ■ , /,»- 2 ) T , f (i) = (/(x.J), • • ■ , /(*S$- 2 )) T , thus 

/j,* = “ 1 < A: < rij — 2. 

i=l 

So we have 

For any given e > 0, we take t = min(rij/ 2, — log e/ log a). For & € [&i + A &2 — f’], using 
(3.28) we have 


(3.28) 


\!.*\ < K{ Y. 3rr ^l/(4? ! )l + E Tiifoil/<*S)U 

a K z l |jk-«l>e aA: z l 


M )| + - 1 
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1 


- XI + ^ ct\k-i\ 

\k-i\<e aK z l |fc-;|>c QA ' Z| 

< 2A'f[l + (-) + h (-)') + 2MK{(-)' + ' + ■ 


a 


= 2A ' f Triir +2M/ ^ ) i-(i) 


< 2/if— ^-r + 2 MKt 1 


a — 1 


a — 1 


= Ce 

where C' = ^^-(a + M). 

Finally, we have 

|L,/(x) - i,j(x) 1 = 1 £ /wnM*)l 

A-etA.-! +/,Ar 2 -0 

< \IjMj,k( x )\ < C'e ^2 \^jA x )\ 

ke[ki+fM-(] kefii+iM-t] 

Note that in the last summation, only three terms will be nonezero for any fixed x, so 
we have 

\I Wj f(x)-l w J{x)\<We = Ce 

where C = + M). This concludes the proof of the Theorem. 


□ 


Remark. As a consequence of Theorem 3, the coefficients fj ^ of the wavelet interpolation 
operator l Wj f(x) can be ignored if x^ € where the function f(x) is less than 

some given error tolerance e. This procedure will only result in an error of 0(c). For 
e = lO -10 ,/ = 9,e = 10 -8 ,/ = 7. In the wavelet interpolation expansion (3.10), I W} is 
used to interpolate the difference between a lower level interpolation Vj-if(x) and f(x), 
i.e. Vj-\J{x) - f(x). Thus, the situation mentioned here will occur in larger region of the 
solution domain as j becomes larger, avoiding adding unnecessary expansion terms ^ h k{x). 
This fact will be used in the later section to achieve adaptivity for the solution of PDE’s. 
The idea of decomposing numerical approximations into different scales has been previously 
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used successfully in the shcok wave computations with uniform high order spectral methods, 
where ENO finite difference methods and spectral methods are combined to resolve the 
shocks and the high frequencey components in the solution, respectively [17]. 

We conclude this section with the following result which shows how to use wavelet coef- 
ficient to estimate the data interpolated by l Wj . 

Theorem 4 Let I W] f(x) and f(x) as in Theorem 4. And if for e > 0, — 1 < k x < k 2 < n 3 —2 

\fj,k\ < f for k x < k < k 2 , 

then 

l/( x **)l < 3c for k x + 3 < k < k 2 — 3. (3.29) 

Proof. The proof follows from the definition of I W} f(x). 

□ 

4 Derivative Matrix V 

The operation of differentiation of functions , which are given in terms of the wavelet ex- 
pansion of (2.39), can be represented by a finite dimension matrix T>. Such matrix has been 
investigated in [16] for wavelet approximation based on Daubechie’s compactly supported 
wavelets for periodic functions. The properties of matrix 2), especially of its eigenvalues, 
affect very much the efficiency and stability of the numerical methods for the solution of 
PDE’s to be discussed in the next section. 

We consider the derivative matrix which approximates the first differential operator 

Cu = u x (4.1) 

with the boundary condition 

u(L) = 0. (4.2) 
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Because of the multiresolution structure of spaces V,-, i.e. Vj C Vj+i and Vo 0 Wo 0 • ■ • 0 
Wj = Vj+i • We can rewrite the wavelet interpolation uj(x ) of (2.39) for function u(x) as a 
linear combination of Ib t j+\f(x) and basis in Vj +\ , namely 

L'- 4 

uj{x) = l h< j+-[u{x) 0 u_i <j>b,j+-i(x) + ^2 iik(f>J+i,k(x) + u L '- 3 <j) bi j + i(L - x) (4.3) 

k= 0 

where V = 2 J+l L and I(,,j + iu(x) is defined in (2.33). 

With the transformation £ = 2‘ /+ 1 x, equation (4.3) becomes 

L '- 4 

uj(£) := uj(x ) = lb,ou(0 + u-i 4>(£) 0 X] Uk<f>k(Q + UL'-3<f>b{L' - 0- ( 4 -4) 

k = o 

Using the notations 


U ' = ( U ( 1 )> U ( 2 ), • * ' yU(L' — 1 )) T G 

a = |.(0],M T ,«(f)) T eR w 
u = (u_i,u 0 , ■ • ■ ,ui-_ 3 ) T e R L ’ - \ 

and equation (3.5), we have 

u = B - 1 (u' - 14 ) 


where vector 14 is defined by 


(4.5) 


u 6 = (I M u( 1 ), 0, • • • , 0, 1 b fi u{L' - 1 )) T € R^'" 1 


and 


l5,0«(l) = ^( c 0 ) 4> C 2> c 3, 0, ■ • • , 0)u ; = 7 iu', 7 , € R 

1 


V - 1 


I 6 , 0 u(L / - 1) = ^(0, ■ • • ,0, -4, -4, -4, -4) = 72U', 72 6 R r i . 


Therefore, 


u = B (I - 


7i 

0 


0 

72 


) U ' = B - 1 ru' 


(4.6) 
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where I is the (Z/ — 1) x (L 1 — 1) identity matrix. 

To obtain approximation to the derivatives of u(£), we differentiate equation (4.4) with 
respect to £ and evalute at fa = k,0 < k < L', i.e. 

*4(6) = M'(W + 

k=0 

= *4(6) + *4(6) 0 <k<V (4.7) 

where u\(fa) denotes the first term in the first equation and u' 2 (fa) the rest. Recalling the 
definition of Ife,ou(£) in (2.33) and coefficients a* in (2.38) with h = 1 and j — J + 1, L = U, 
we have 


*4(0) 

*40) 


2ELo c ’k u i k ) - MO) 
-m=oc'Mk) 


x lES-o c'ML'-k) 

\ -2n= 0 c'MU-k) + 3u(L') 

with the four L' + 1 dimension vectors 



u = Au 


= (24 - 3, 2 C ; , 2c', 2c' , 0, • • • , 0) € R L ' +1 
= -^(4,c' 1 ,c',c',0,---,0)eR i ' +1 
= ^(0,---,0,C3,C2,c;,Co) € R L ' +1 
= -(0, • • • , 0, 2c', 2c' , 24,24 - 3) 6 R L ' +1 . 


On the other hand, using (3.5) we have 


0 0 


*4(0) 

*40) 

*4(M 


000 


0 0 
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(4.9) 


= HB-'Tu' 

= (o,^B- 1 r,o)u 

Finally, combining equations (4.8) and (4.9), we have 

/ D^uj( 0 ) 

^TIj(l) 
u < = : 

\ D ( uj{L') 

where the derivative matrix T>' is defined by 

P' = A + (0,//B" 1 r,0). (4.11) 

Converting to the x- derivatives, we have 

D x uj(x 0 ) \ / D^uj( 0) 

D x uj(x i) _ J+1 D ( 'Uj{ 1) 

D x uj{x L ') / \ D^uj(L') 

where x; = :pTr> 0 < i < Z/. 

Let P be the upper left L' x L' submatrix of P', will be the wavelet derivative 

matrix to differential operator (4.1) with boundary condition (4.2), namely, P will maps the 
function values u(xo)) u(xi), ■ ■ • , u(x£/_i ) to its derivatives u^Xo), u'(xi), ■ ■ ■ , u (xz,'_i), 

u'(x 0 ) \ / u(x 0 ) 

“'(*>) =2 J+, P “ (l,) 

u'(xl>- i) / \ u(xl'-i) 

In Figure 15, we plot the eigenvalues of P for L = 8, J = 0 , 1,2,3 which corresponds to 

N = 8,16,32,64. The eigenvalues come in conjugate pairs with two pure real eigenvalues. 
The real part of all the eigenvalues are negative and except one eigenvalues, all the rest are 
close the imaginary axis. 

5 Adaptive Wavelet Collocation Methods for PDE’s 

In this section we consider a collocation method based on the DWT transform given in 
Section 3 for time dependent PDE’s. Let u = u(x,t) be the solution of the following initial 
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value problem 


«t + fx(u) = u xx + g(u), x e [0, L\, t > 0 
u(0,t) = g 0 (t) 

u(L,t) = 9l (t) l° A > 

< u(x,0) = f(x) 

where only Direchlet boundary conditions are considered, however, the methods presented 
here can also be modified to treat Von Neuman type or Robin type boundary conditions. 

We use the idea of method of lines where only the spatial derivative is discretized by the 
wavelet decomposition. The numerical solution uj(x,t ) will be represented by an unique 
decomposition in Vo ® W 0 0 • • • © Wj, J > 0, namely 

L— 4 

uj(x,t) = \ju(x,t) + U_l -i(t)(f> b (x) + + U^ itL _ 3 (t)(f> b (L - x) 

k=0 

, 7=0 k=— 1 

J 

= u-i(x) + (5.2) 

3=0 

where I bl ju(x,t) given in (2.33) consists the nonhomogenuity of u(x,f) on both boundaries, 
and the coefficients Uj^t) are all functions of t. Using the DWT transform, we can also 
identify the numerical solution uj(x,t ) by its point values on all collocation (previously 
named interpolation ) points, i.e. {arj^} in (3.1) and (3.6), we put all these values in vector 
u = u (/), i.e. 

u = u(i) = (u^ -1 \ u*°\ • ■ - ,u^) T 

where = {u(xjj. j) , f)}, 1 < k < L - 1 for j = -1; -1 < k < n, - 2, for j > 0. 

To solve for the unknown solution vector u(/), we collocate the PDE (5.1) on all colloca- 
tion points, then we have the following semi-discretized wavelet collocation method. 

Semi-Discretized Wavelet Collocation Methods 

ujt + fx(uj) = Uj xx + g(uj)\ x=x (j),—\ < j < J 
j U A <M) = 3o(t) .p 

uj(L,t ) =</i(t) 

„ u j{ x = x[ j) ,0) = f(x = x^) 

where 1 < k < L — 1 for j = — 1 ; — 1 < k < tij — 2, for j > 0 


24 


iiUiin iin iluiimi ii iinii n ii m m 



Equation (5.3) involves total (2 J+1 - \)L + 2 unknowns in u two of which will be de- 
termined by the boundary conditions and the rest are the solutions of the ODE system 
subject to their initial conditions. In order to implement the time marching scheme for the 
ODE’s system (for example Runge-Kutta type time integrator), we have to compute the 
derivative term in (5.3) f x (uj(x[^)) and u Jxx (x [ 3) ) in an efficient way. Let us only discuss 
the first derivative which involves the computation of the nonlinear function f(uj(x,t)). For 
this purpose we first find a similar wavelet decomposition as (5.2) for f(uj). For a general 
nonlinear function /(u), this can be done quite straightforward using the DWT transform 
in Section 3. 

Computation of / X ((xj^) = f x (uj(x[ p) 

Step 1 Given u = (u (-1 \u (0) , ■ • • ,u (i) ) T , compute f (j) = j > -1 and define 

f = (fM) ) f(°) ? ... ) fG ) ) T ; 


Step 2 Compute the wavelet interpolation expansion using DWT transform for f, 


L - 4 
k=0 

+ Et E *(*)]; 

j= 0 fc=— 1 

Step 3 Differentiate (5.4) and evaluate at all collocation points > —1, 


(5.4) 


Muj) i x _ x (» = (i m/)'( 4 j) ) + /- i.-i + E f-ukW k (4 j) ) - -4 J 

* Jt=0 

+ Et E 

i=0 (=-l 


U) 


Cost of Computing the Derivatives. 

For each single collocation point, it takes 7 + 5 (J + 1) = 5 J + 12(//op.s) to compute 
Therefore, the total cost of computing all derivatives is (5J + 12) N < 5NlogN. 


Again, t/>'(x) and <f>'(x) at the dydic. points 0 < k < 2* L can be precomputed once and 
for all. 

Assuming that Euler forward is used to discretize the time derivative in (5.3), we obtain 
a fully discretized wavelet collocation method. 

Fully discretized Wavelet Collocation Method 

( uj n+ ' = uj n + At[-f x (u n j) + u n Jxx + g(u n j)}\ x = x oh -1 < j < J 

< «3(0) =90 (t n ) 

u j(L) =9*(t n ) 

. u °j( x = 4 J) ) = f( x = x i 3) ) 

where 1 < k < L — 1 for j = — 1; — 1 < k < nj — 2, for j > 0 and t n = nAt is the time 

station and At is the time step. 

Adaptive Choice of Collocation Points 

* 

In equations (5.2) and (5.4), uj(x ) and f(uj(x)) are expressed using the full set of 
collocation points {#* }. As discussed in the remark after Theorem 3 of Section 3, most 
of the wavelet expansion coefficients uy* for large j can be ignored within a given tolerance 
e. So we can dynamically adjust the number and locations of the collocation points used in 
the wavelet expansions, thus reducing significantly the cost of the scheme while providing 
enough resolution in the regions where solution varies much. We can achieve this adaptivity 
in the following two ways. 

Deleting Collocation Points 

Let f > 0 be a prescribed tolerance and j > 0, i = £(e) = min ( — log e/ log a). 

Step 1. First we locate the range for the index k, 

(k'j,lj),---,(k^,l' m ),m = (5.6) 

such that 

< (, K < k < i = 1, • • • ,m. (5.7) 


(5.5) 
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Step 2. Following Theorem 3 and 4, we can ignore iijj t in u 3 (x) in (5.2) for k, < k < 
/, , i = 1 , ■ - ■ , m, ki = fc t - -f ^ + 3 ^ — 3, namely we redefine u 3 (x) as 


where Kj = Ui <.•<«[*■, *,]• 

Step 3. The new collocation points and unknowns will be 

= 1, 1 if j = -l;fc € 1, — - 2}\/Cj. 


Increasing Level of Wavelet Space. 

Let t > 0 again be some prescribed tolerance, and if 


max\u Jk \ > e 


(5.8) 


where subscript n indicating the solution at time t = then we can increase the number 
of wavelet spaces W 3 in the expansion for the numerical solution uj(x) in (5.2), say, up to 

Step 1 At t = t n if condition (5.8) is satisfied, let J' > J and define a new solution vector 

u 3, := (u (-1) , U W, ■ • ■ , uW, • • • , u<- , ')) T 

where for J + 1 < j < J\ = {'Pju(x^)}^~ 2 1 . 

Step 2 Use u”, on the right hand side of scheme (5.5) to advance the solution to time 
step f n+1 and obtain solution Uji*’ 1 . Then, Ujt*(x) = Vj>Ujt 1 6 Vo © Ho ® • • • ® Wj> will be 
the new numerical solution which yields better approximation to the exact solution of (5.1). 

6 Numerical Results 


CPU Performance of DWT transform 

The theoretical estimates of operations for performing the DWT transform in both di- 
rection and the computation of derivatives at all collocation points are 0( N log N ) where N 
is the total number of terms in the wavelet expansion (3.10). 
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We take the function in (6.1) and define its wavelet interpolation expansion (3.10) for 
L = 10, J = 2, 3, • • • , 9, the total number of terms (or collocation points) N = 2 7+1 L - 1 
are between 79 and 10240. In Figure 5 we plot the CPU time for the performance of DWT 
back and forth in both directions (‘o’ in the Figure) and the computations of derivatives on 
all collocation points ('+’ in the Figure). Also drawn in the Figure is a straight line which 
indicates a almost linear growth of the CPU timing up to 10 k points. 

Adaptive Approximation of Wavelet Interpolation Expansion 

We consider function- 


/(*) = 


[ h^(x + 1,0.3) 

0 

h\{x + 0.5, 6) 

0 


if - 1 < x < -0.7 
if -0.7 < x < —0.5 — 8 
if - 0.5 - 8 < x < -0.5 + 6 
if —0.5 + 6 < x < 0 


( 6 . 1 ) 


sin(57rx)A,(x - 0.25,0,25) if 0 < x < 0.5 
{ h 2 (^) if 0.5 < x < 1 


where 6 = 0.01 and hi(x,a ) is an exponential hat function and h 2 {x) is a step-like function 
and they are defined as 


-J 

exp( a2 l x2 ) 

if |x| < a 

(6.2) 

l 

0 

otherwise 

0 


if x < 0 


1 

2772 

J, o* 5 (l-<) 5 

dt if o < x < i 

(6.3) 

1 


otherwise 



and 


First we construct the full wavelet interpolation expansion (3.10) Vjf{x) for J = 6, L = 
40, the total number of wavelet functions (or the collocation points N ) N + 4 = (2 J+l L - 
1) + 4 = 2 J+1 L + 3 = 5123 (including four boundary functions in 1 6, •//(*))• Figure 6, on 

the top we plot the f(x) (solid line) and Vjf{x) at non-interpolation points, at the bottom 
we have the absolute error in logarithm scale. In Figure 7, we plot the components /o G Vo 
and 9 j (x) € Wj, 0 < j < 6 in Vjf(x) = I b ,jf{x) + /o + 9o + • • • + 9J- We can see that only 
higher frequency part is retained in higher wavelet spaces Wj (notice that the scales varies 
in different pictures). 
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Then, we use the procedure at the end of Section 5 to filter out those coefficients, thus 
deleting the corresponding collocation points, fj t k which are less than e in magnitude. In 
Figure 8, we take c = 10~ 5 and the number of wavelet functions reduced to 289 with the 
accuracy of the approximation (bottom curves) within order of e. In Figure 9, we plot the 
solution at the remaining interpolation points and the expected clustering of the interpolation 
points is seen at locations where the function changes more dramatically. In Figure 10, we 
plot the magnitude of the wavelet coefficients fj f k,j > —1 one level above another. High 
density of the wavelet coefficients reflects the existence of high gradients of the approximated 
function. In Figure 11, we take e = 10 -4 and the number of wavelet functions fj t k reduced 
to 206 with the accuracy of the approximation (bottom curves) within order of e. 

Linear Hyperbolic PDE’s 

We consider the IVB problem of linear hyperbolic partial differential equation 

ut + u x = 0, 0 < x < 1 

< u(0,0 = 0 (6.4) 

. «(z>0) = h 2 (§) 

where 8 = 0.05 and h 2 {x) is defined in (6.3). 

We apply the collocation method with adaptive choice of the collocation points L = 
20, J = 4. Second order Runge-Kutta method is usd for the time derivative. With every 
10 iterations we change the number and locations of the collocation points according to the 
criteria proposed at the end of Section 5. The cut-off tolerance e = 10 -5 . The number of 
collocation points involved fluctuates around 200 in contrast to the full set collocation count 
which is 640 in this case. In Figure 12, we plot the numerical solution (' + ’) against the exact 
solution (‘o’) at time t = 0.1. In Figure 13, we plot the errors in logarithm scale (notice the 
y-sc.ale starts at -2 which corresponds to an error of 10~ 2 ). Again, we see the automatically 
clustering of the collocation points. 

Inviscid Burger Equation 
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Finally, we consider the IVB problem of nonlinear hyperbolic partial differential equation 

f «« + (£)* = <>, -1 <* <2 

< u(0,<) = given (®-®) 

[ u(x, 0) = f(x) 

where 

f lsin{-Kx) if — 1 < x < 1 
J\ x ) ~ | q otherwise 

In this case, we take L = 10, J — 6. Second order Runge-Kutta method is used for 
the time derivative. With every 10 iterations we change the number and locations of the. 
collocation points according to the criteria proposed at the end of Section <). The number 
of collocation points involved fluctuates around 100 in contrast to the full set collocation 
count which is 1280 in this case. The cut-off tolerance e = 10~ 5 . In Figure 14, we plot the 
numerical solutions at time t = 0.05,0.1. The numerical scheme automatically puts more 
collocation points near the high gradient (x=0) and the derivative discontinuity (x=l). 

7 Conclusion 

In this paper, we have constructed a fast Discrete Wavelet Transform (DWT) which enables 
us to study collocation methods for nonlinear PDE’s. The adaptivity of wavelet approxima- 
tion is conveniently implemented through the examination of the wavelet coefficients. The 
preliminary tests on the solution of PDE's indicates such an approach will be important 
in large scale computation where the solution develops extremely high gradients in isolated 
regions, and uniform mesh is not practical. Such investigations are actually being done for 
reacting flows, the results will be reported in a separate paper. 
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Appendix 


Proof of (3.25). The proof is a straightforward application of Lemma 5. For M } in 
(3.9), we have (a 2 , a 3 , • • • , a n ) = (-i • ■ • , -£), (6,, & 2 , • • ■ , fe n ) = (1, 1, ■ ■ • , 1) and 

( c 2, c,3, • • • , c,,) = (— P4, — pj, ■ ■ • , — yj,— pj) where n — nj = 2 PL. Therefore, the sequence 
{u TO } in ( 3 . 22 ) satisfies the following relations, 


«o = 0, = 1, u 2 = 14, u 3 = 


2534 
13 ’ 


(A.l) 


and for 4 < m < n 

Mm = Cm{ 2 4" 14u m _j) (A. 2) 

where c m = if m = n, c TO = 1 otherwise. 

Recursive relation (A. 2) is a finite difference of order 2 whose general solution is of the 
following form 

«rn = c, n (c,a m-3 + c 2 /? m - 3 ) (A. 3) 

where cv = 7 + \/l 92/2, (3 = 1 — 1 92/2 are the two distinct roots of the quadratic equation 


2? 2 — 14x T 1 — 0, 


and constant cj and c 2 are chosen so equation (A. 3) is valid for m = 2,3. 
Therefore, 


u 0 = 0,i/i = 1 

Mm = c m (//,a m “ 3 + / / 2 /? m " 3 ), 2 < m < n (A.4) 

where //, = ^(^ - 14/3) > 0,// 2 = ^(14ar - 2524) > 0. 

Similarly, we can show that 

Mn+l = 0, Vn= 1, (A. 5) 

v m — c n _ m+ i(//]O n ~ 2-TO -f- fi 2 /3 n ~ 2 ~ m ), for 1 < m < n — 1 (A. 6) 
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and 


»o = -!(«, a"" 1 + V*). 

(I 1 

where <5, = < 13a ~ - ^ > 0, <5 2 = < 0 . 

Finally, following (3.24) in Lemma 5, we have the following estimates on the 

Mj. 

Denote ej, 1 < j < n as 

r i if j = i 


e i = 


8 if ^=2 


I if 3 < j < n - 1 

II 
13 


li if j = n. 


Case 1 : i < j and 1 < j < n — 1 , i = 1 


&ij — ^i^n— j+\ 


{6, a n ~ 4 + 6 2 /? n ~ 4 ) 


So we have 


l®ijl < 


< 


14 a n ~ 2 ~ J (/ui +A t 2^ n ~ 2 ~ J ) 

13 a n_4 (5i — S2Z n ~ 4 ) 

14a(/^i + / g ) 1 


13 (<$i — 2) or J 1 


Qfl J — *1 


where z = ^ < 1 and A^i — 1.1666. 

Of — 

Case 2: i < j and 1 < j <n — 1 , 2 < i < n 


_ (//ia ; - 3 + M2/? , ~ 3 )(^i« n ~ 2 ~ i + M2^ w - 2 ~ j ) 
e * c "-i+ lC ‘ (^a 1 *- 4 +£ 2 /? n " 4 ) 


Case 3: i < j and j = n,i = 1 


1 

a i,J ~ -e, (^ a n-4 + 8 2 {l n - 4 ) 


Case 4: i < j and j = n, 2 < i < n 

- (f 1 l a ' _3 + ^2^' 3 ) 
aiJ = ' e, ' Ci (« 1 a"- 4 +«j|9"- 4 )' 


inverse of 


(A.7) 


(A. 8) 


(A.9) 


(A. 10) 
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■•I 'll ■inum mil M * hi i mm nikN p || || , I. I,I. I. II , . M nil 



Case 5: i > j and j — 1,1 < i < n — 1 


a 


hj 


= —tjC 


«-t+l 


($ia n - 4 + M n ” 4 ) 


Case 6: i > j and j = 1 , i = 


71 


~ e j 


(8ia n ~ 4 + 6 2 (3 n ~ 4 ) 


Case 7: i > j and 2 < j < n — 1, 1 < i < n — 1 

(^i<* J ~ 3 + ^ 3 ~ 3 )(fi l a n - 2 ~ t + fi 2 /3 n ~ 2 ~') 


C jCjC n —{_ j_] 


Case 8: i > j and 2 < j < n — 1 , i = n 


a i,j — e 


(V’ i - 4 + ^ n - 4 ) 


(Hia j 3 + fi 2 (3 j 3 ) 


For Cases 2 - 8, we can similarly obtain 

K, 




a 


li-'l ' 


2 <i <8 


where K 2 = 1.1726, K z = 1.1607, = 1.1666, A's = 1.1666, /C 6 = 1.1607, K 7 = 

h' s = 1.1666. 

Finally, if we choose K = 1.1726, then 

. « K 

Kil £ -\ ifTr 1 ^ < »• 


This concludes the proof. 


(A. 1 1 ) 


(A- 12) 


(A- 13) 


(A- 14) 


1.1722 and 


(A. 15) 
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Figure l.Interior scaling functions 4>{x) (top) and boundary scaling function 
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Figure 2. Interior wavelet functions rj>(x) (top) and boundary wavelet function if> b (x). 
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Figure 3. Fourier Transformations of x(){x) 
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Figure 4. Fourier Transformations of ip"(x) 
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Figure 5. CPU timing for Performing DWT (both directions) transformation (‘o’) and 
computation of derivatives (*+’), solid line - Linear fitting. 
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Figure 6. Wavelet approximation of function (6.1) with L = 40, J = 6. Top - Exact 
solution (solid line) and approximation (‘o’); Bottom - absolute error in logarithm scale. 
Total number of is 5123. 
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Figure 8. Same as Figure 6, but with deletion of wavelet coefficient whose magnitude 
less that e = 10 -s . Total number of fj^ left is 289. 
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1 < 

0.8 

0.6 

0.4 

0.2 



Figure 9. Close up of top part of Figure 8, numerical solutions (*+’) at remaining collo- 
cation points against exact solutions ('o’). 
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Figure 10. The magnitude of remaining wavelet coefficient in Figure 8. 
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- 0.8 - 0.6 - 0.4 - 0.2 0 0.2 0.4 0.6 0.8 


Figure 11. Same as Figure 8, but with e = 10 -s . Total number of fj,k left is 206. 
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-1 - 0.5 0 0.5 1 


Figure 12. Adaptive collocation solution of linear PDE (6.4) at t = 0.1 with L = 20, J = 4 
and error tolerance e = 10 -4 . Total number of collocation points is around 200. 
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^ ■ ■ > , i . 7 

-1 - 0.5 0 0.5 1 1.5 2 


Figure 14. Adaptive collocation solution of non-linear PDE (6.5) at t = 0.05, 0.1 with 
= 10, J — 6 and error tolerance e = 10~ 5 . Total number of collocation points is around 
100. 
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Figure 15 Eigenvalues for the first derivative V for L = 8,7 = 0, 1,2,3 whose sizes are 
2 J L — 8,16,32, and 64, respectively. 
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